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ABSTRACT 

We analyse the interaction of an eccentric binary with a circular coplanar circumbinary disc 
that rotates in a retrograde sense with respect to the binary. In the circular binary case, no 
Lindblad resonances lie within the disc and no Lindblad resonant torques are produced, as 
was previously known. By analytic means, we show that when the binary orbit is eccentric, 
there exist components of the gravitational potential of the binary which rotate in a retrograde 
sense to the binary orbit and so rotate progradely with respect to this disc, allowing a resonant 
interaction to occur between the binary and the disc. The resulting resonant torques distinctly 
alter the disc response from the circular binary case. We describe results of three-dimensional 
hydrodynamic simulations to explore this effect and categorise the response of the disc in 
terms of modes whose strengths vary as a function of binary mass ratio and eccentricity. These 
mode strengths are weak compared to the largest mode strengths expected in the prograde 
case where the binary and disc rotate in the same sense. However, for sufficiently high binary 
eccentricity, resonant torques open a gap in a retrograde circumbinary disc, while permitting 
gas inflow on to the binary via gas streams. The inflow results in a time varying accretion rate 
on to the binary that is modulated over the binary orbital period, as was previously found to 
occur in the prograde case. 
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1 INTRODUCTION 

We investigate the behaviour of retrograde circumbinary discs for eccentric orbit binaries. When the binary orbit is circular, retrograde 
circumbinary discs can behave fundamentally differently to their prograde counterparts, due to the absence of orbital resonances (Nixon et al. 
2011a). In the prograde binary case, resonances between the binary and disc transfer angular momentum from the binary to the disc. This 
transfer causes the binary orbit to shrink somewhat, perhaps also changing the binary eccentricity, while the disc gains angular momentum 
and experiences gap opening (see the review by Lubow & Artymowicz 2000). In a one dimensional idealisation, the prograde disc is therefore 
often described as a decretion disc (Pringle 1991) in which there is a central torque and no gas flow onto the central binary, rather than a 
standard accretion disc in which there is no central torque and there is gas flow onto the central object (Pringle & Rees 1972; Pringle 1981). 
However, this idealisation of a decretion disc is often not realized in multi-dimensional studies. Although a tidally produced gap does form, 
gas flows through the gap in the form of streams, resulting in substantial accretion onto the binary (Artymowicz & Lubow 1996; Gunther & 
Kley 2002; MacFadyen & Milosavljevic 2008; Shi et al. 2012). The accretion rates can be comparable to the rates expected from an accretion 
disc that surrounds a point mass whose mass is equal to the binary mass. The gap is formed by the resonant interaction of the binary with the 
disc (Artymowicz & Lubow 1994). 

In contrast to the prograde case, a retrograde circumbinary disc experiences no resonances and therefore flows towards the binary as 
a standard accretion disc until the disc orbits are perturbed significantly by close passages with the binary. The perturbed material is then 
captured into circumprimary or circumsecondary discs (Nixon et al. 2011a). For a circular orbit binary in a coplanar circumbinary Keplerian 
disc, Lindblad resonances occur where 

(r) = m\a{r) - , (1) 
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where Hb is the binary orbital frequency, r is the distance from the binary center of mass in the binary orbit plane, Q(r) is the disc orbital 
frequency, and m > 0 is the integer wave mode number (e.g., Goldreich & Tremaine 1979). For a retrograde disc, the two frequencies have 
opposite signs and the right hand side of Equation (1) is always larger than the left hand side. Therefore there can be no resonances in the 
disc. 

In the case of eccentric orbit binaries, many more resonances are produced whose associated torques depend on the binary eccentricity 
(Goldreich & Tremaine 1980). Such resonances are capable of opening central gaps in prograde coplanar circumbinary discs (Artymowicz & 
Lubow 1994), while allowing gas flow through the gap in the form of time-modulated gas streams (Artymowicz & Lubow 1996). However, 
as we shall see below, when the binary is eccentric, it is possible to generate resonances between the binary and a counter-rotating disc. This 
effect is due to components of the eccentric binary potential that rotate in the opposite sense to the binary, and therefore progradely with 
respect to the disc. We describe this effect analytically and explore it with three dimensional hydrodynamic simulations. 

Retrograde circumbinary discs (Nixon et al. 2011a; Nixon 2012; Roedig & Sesana 2014) have received relatively little attention com¬ 
pared to prograde discs (e.g. Artymowicz & Lubow 1996; Escala et al. 2005; Lodato et al. 2009; Cuadra et al. 2009; Roedig et al. 2012; 
D’Orazio et al. 2013). The simulations of Nixon et al. (2012) and Nixon (2012) focussed on circular orbit binaries, with only a few sim¬ 
ulations with ^ ^ 0, in which the eccentricity was always small. Roedig & Sesana (2014) performed simulations of retrograde discs with 
eccentricities as high as ^ = 0.9, but due to the strongly self-gravitating discs employed, and the induction of disc tilt due to the dominant 
disc angular momentum (cf. King et al. 2005; Nixon et al. 2011b), the subtle resonant effects sought here may have been missed. In this 
paper we focus on the case where the disc mass is negligible and so the binary dominates the system angular momentum, even in the highest 
eccentricity cases. 

We present an analytic model and simulations of retrograde circumbinary discs, varying both the binary mass ratio and eccentricity. In 
Section 2 we define mode strengths that are an extension of the definition given by Lubow (1991) to handle retrograde discs. In Section 3, 
we evaluate by means of linear theory the torques produced in a retrograde disc due to an eccentric orbit binary. Section 4 describes a set of 
simulations and their agreement with the analytic model. Section 5 contains a discussion, and Section 6 contains the conclusions. 


2 MODE STRENGTHS 


In order to gain insight into the response of a disc to tidal forcing, we decompose this response into different modes that are characterized by 
azimuthal mode numbers m and frequency mode numbers We consider a surface density distribution E(r, 0, t) and write 

oo oo 

S(r, e,t)=Y,Yj exp [i{m0 - (t)] ], ( 2 ) 

i=-oo m=0 

where E^ ,„(r) is a complex function and r = is the mean anomaly of the binary orbit that may be eccentric, with mean motion Ob > 0 
(defined as 2;r/Pb for binary orbital period Pb)- We invert Equation (2) to obtain 

^2n ^2n 

'Le„(r)=—^ ——r I I S(r, 0,0 exp [-((mS - fr)] Jt. ( 3 ) 

2;r^(l -I- O{fidm,o) Jo Jo 

Writing the integral in Equation (3) in terms of real functions, we obtain 

^ejn(r) = ^cos,cos/,m(r) + ^sin,sin/,m(r) + (^sin,cos/,m(r) ^cos,sin/,m(r)) (4) 


^ ^2n ^2n 

c c ^ I I ^(r,e,t)f(m 0 )gi£T)d 0 dT, 

^ O£,0Om,0) Jo Jo 


where 

2f,g/,mW= 

and / and g each can be the cos and sin functions. We define a radially integrated quantity X as X so that 

^oo 

— I f'dr — 2icos,cos/,m ^sin,sin/,m + (^sin,cos/,m “ ^cos,sin/,m)9 

Jo 

where 




0 




We then determine the square of the radially integrated mode amplitudes that are given by as 


~ (^cos,cos/,m ^sin,sin/,m) (^sin,cos/,m ^cos,sin/,m) 


T + 

cos,cos/,m sin,sin,^,m cos,sin,^,m sin,cos,^,m 


2 ^^cos,cos,^,m^sin,sin,/’,m ^sin,cos,^,m^cos,sin,/’,m^ • 


> e,m ■ 


We define the dimensionless mode strength S such that S o,o is unity 

Mi ’ 


(5) 

(6) 

(7) 

( 8 ) 

(9) 
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where is the disc mass. This definition of mode strength is similar to that given in Lubow (1991). One difference is that the definition 
in Lubow (1991) included only the squared terms in the second line of Equation ( 8 ) and omitted the cross terms. To see the effect of the 
omission, consider the case of density distribution 

2:(r, 6 >, t) = S(r - 1) cos (m'6 - rObO- (10) 


In that case, we have that 

(H) 


Using only the squared terms on the second line of Equation ( 8 ), as in Lubow (1991), we find that 

~ 2 ^Scos,cos/,m^sin,sin/,m ~ ^sin,cos/,m^cos,sin/,m^ “ • 


( 12 ) 


Equation (12) shows that the previous definition of mode strengths does not distinguish between prograde and retrograde modes in which 
the sign of £ differs. It was adequate for the analysis of the superhump instability cycle that only involved prograde (£ > 0) modes, but is 
not adequate for the current study. Instead we must use the definition of given by Equation (9). In Appendix B, we describe how mode 
strengths are computed in SPH simulations. 


3 RESONANT TORQUES 


3.1 Torque equation 


We determine the resonant torque on a gas disc due to an eccentric binary that is exerted on a counter-rotating disc. The disc mass is assumed 
to be very small compared to the binary mass. We take the binary to be on a prograde orbit with mean motion Ob > 0, while the disc orbit is 
prograde 0(r) > 0 or retrograde Q(r) < 0 with respect to the binary. We adopt a cylindrical coordinate system whose origin is at the binary 
center of mass and so the z = 0 plane coincides with the binary orbit plane. Following Goldreich & Tremaine (1979), we decompose the 
potential as 

0(r, 6,t) cos {mO - ^QbO, (13) 

i,m 

where m > 0 and i ranges over all integers (negative, zero, and positive). The binary potential has contributions from the primary star of mass 
Ml and the secondary star of mass M 2 . We determine the potential components for £ and m > 0, by 

^eAr) = + <D® (r) (14) 


for m > 0 where the superscripts 1 and 2 denote the primary and secondary objects. Following Moorhead & Adams (2008), we invert 
Equation (13) by using the eccentric anomaly of the binary orbit ^ as a variable of integration in place of the mean anomaly, Q\yt. We then 
have for / = 1,2 


O 


(i) 

i,m 


(r) = - 


GMi 



(1 - e cos ^) cos (mO - £(( - e sin ^)) 
+ xf (I - e cos - ipxtgiO, ^) 


(15) 


where 


^) = (cos ^ - e) cos {6) + Vl - sin (f) sin {6), 


(16) 


xi - -M 2 1 {Ml -h M 2 ) and ^2 = Mi/(Mi -h M 2 ), e is the binary eccentricity, and p - rja, for a binary with semi-major axis a. The advantage 
of this method is that we can determine numerically for arbitrary binary eccentricity. We have verified analytically in expansions for 
small e that this method recovers potentials Oi 1 , Oi 2 , and Oi 3 given by Equations (18), (21), and (23) respectively of Artymowicz & Lubow 
(1994). 

We determine the resonant torque on the circumbinary disc due to outer Lindblad resonances that occur where 


mfl(r) - ffib = -K{r), 


(17) 


where k is the epicyclic frequency. Note that we do not assume the disc is Keplerian. We take into account the non-Keplerian effects of the 
binary gravitational potential by determining Q(r) and /c(r) as 


a\r) = l^ 

r dr 


(18) 


and 


2 3 dOofi d^Oo^o 

where we take Q and k to be negative for a counter-rotating disc. 
In the case of a Keplerian disc, we have that CI-k and so 


( 19 ) 
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n(r) = 




( 20 ) 


m + 1 

at an outer Lindblad resonance. This estimate treats the binary as a point mass which is valid at large distances from binary. However, as 
noted above, we do not make this approximation in the torque calculations that we carry out. 

The magnitude of the resonant torque on the disc (Goldreich & Tremaine 1979) is given by 




\D\ ’ 


where 


= r- 




c,m 


2mCl 


-o, 


dr mQ(r) - ffib 


i,m 


and 


D = r 


d{K^ - (ma - 
dr 


( 21 ) 


( 22 ) 


(23) 


In Equations (21) - (23), all radially dependent quantities are evaluated at the resonance radius given by Equation (17). We find that the outer 
Lindblad resonance torques we consider lead to pushing material outward for both prograde and retrograde discs. 


3.2 Results 

We describe here the results of the application of the torque model to the cases that we have analysed in the simulations discussed later in 
Section 4 that have values of (1,1), (1,2), (-1,0), (-1,1), and (-1,2). Erom Equation (20), it is then clear that torques with negative £ 
values are possible in counter-rotating discs for eccentric binaries, since Q is negative at resonance. However, since m > 0, it then follows 
that such resonances involve £< 0 < m. Goldreich & Tremaine (1980) showed that 7/,^ oc for ^ 1. That is, resonant torques are 

possible in a counter-rotating disc for eccentric binaries (e > 0 ) provided that ^ < 0 < m, but they are generally weaker for larger differences 
between m and £. We then expect that the strongest Lindblad torques that are produced in a disc that counter-rotates with respect to the sense 
of orbital rotation of a binary to have £ = -I. Eor this reason, we have concentrated on resonances with £ = -1, 

Eor the cases of disc modes (£, m) of (1,1) and (1,2), we see from the approximate Equation (20) that resonances occur only for G(r) > 0. 
Since the simulated discs are retrograde Q < 0, we expect that there are no resonant torques for these simulated cases. Therefore, we predict 
that S 14 and S 14 are both zero. These cases serve as as useful check on the both the theory and the numerical accuracy of the mode strength 
determinations from simulations. 

Eor the case of disc mode (-1,1), Equation (20) estimates that the resonance is located at r = 1.58a. This resonance is sufficiently close 
to the binary that nonKeplerian effects can be important. In addition, the binary may cause the disc to be dynamically unstable. Consider first 
the case of a circular binary. If we approximate the disc as consisting of simple periodic orbits in the corotating frame of the binary, we can 
apply the previously known results on ballistic particle orbits in the restricted three-body problem to study the existence and stability of disc 
streamlines near the resonance. A similar approach for the prograde case with a circular binary was taken by Paczynski (1977). Circumbinary, 
retrograde, simple periodic orbits are known to be quite stable, even if they come close to the binary, although the orbits become increasingly 
noncircular close to the binary (Szebehely 1967). We have found that there are only mild departures from circular orbits near this resonance 
location in the case of a circular binary. 

In the case that the binary is eccentric, the situation is somewhat different. Due to the non-Keplerian effects of the binary, the epicyclic 
frequency /c is a function of binary eccentricity. At locations close to the binary, k can vanish for sufficiently high binary eccentricity. In Eigure 
1 , we plot as heavy lines the critical values of binary eccentricity ^cr for which k vanishes as a function of distance from the binary center 
of mass. The plot shows that is smaller closer to the binary, as one would expect. Eor eccentricities above the line, is negative which 
means that the epicyclic motions are unstable and resonance is not possible. The equal mass binary case is more stable than the unequal mass 
cases, possibly because the binary members are located farther inward from a given radius in the circumbinary disc. 

Also in Eigure 1, we plot in light lines the critical eccentricity for which the apocentre of the secondary object reaches a given radius. 
Above these lines, the secondary member of the binary would cross the given radius. Consequently, disc-like orbits should not exist at such 
radii. Eor an equal mass binary, the apocentre reaches a radius rja = I when e = \. The binary cannot extend beyond rla> I for any value 
of e < 1. Consequently, we do not plot a light solid line in the figure. The light lines lie above the corresponding heavy lines. The plot then 
shows that for such critical binary eccentricities with 1 > ^ > 0 . 1 , the disc orbit is already unstable, since < 0 . 

The angular speed G(r), as well as /c(r), is affected by the binary eccentricity. As a consequence, the resonance location given by 
Equation (17) changes with binary eccentricity. Eor the (-1,1) resonance, we find that the effect of binary eccentricity is to pull the resonance 
inward of the Keplerian location of 1.58a. As a result, the critical value for the eccentricity at this resonance is even smaller than would be 
predicted by Eigure 1 for the Keplerian location. We find that the critical value of eccentricity for a binary of mass ratio q = 0.3 is ^cr - 0.55. 
Above the value, we expect r_i 1 and S-i^i to vanish for q = 0.3. In addition, from symmetry considerations, it it easy to see that 0 _i 4 = 0 
for ^ = 1. Therefore, T-ij and 5 -14 should vanish in the case of an equal mass binary. 

Similar stability considerations apply to the resonance for (-1,0) that lies even closer to the binary. In this case, the torque is always zero 
because m = 0 for this mode. However, an energy flux can be generated at such a resonance and S'- 1,0 can in principle be nonzero. Equation 
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Figure 1. The heavy lines plot critical eccentricity values for which the epicyclic frequency k vanishes as a function of radius. The light lines plot the critical 
eccentricity values for which the apocentre of the secondary object reaches a radius on the horizontal axis. The solid, dashed, and dotted lines are for binary 
mass ratios of ^ = I, q = 0.3, and ^ = 0.1, respectively. No light solid line is plotted because a member of an equal mass binary cannot reach the plotted radii 
(r/fl > 1) for any value of eccentricity. 


(20) suggests that the resonance occurs where H = -Qb. if the disc were Keplerian there. Due to the strong nonKeplerian effects, we expect 
that the resonance condition is even more difficult to satisfy in this case than in the (-1,1) case. Consequently, we expect S to also be 
small. 

Torque r_i ,2 occurs near Q(r) = -Hb/S or r ^ 2.08a that is far enough from the binary to permit resonances to operate. At this 
resonance, the disc is in nearly Keplerian motion. In this case, for small e, we expect that r_p 2 ^ oc and the torque is consequently 
very sensitive to the binary eccentricity. In a prograde disc, torque 7^2 occurs at the same radial location. The prograde counterpart to torque 
r_i 2 in a retrograde disc is then Ti^. This torque ri 2 is frequently found to truncate the inner regions of prograde discs that orbit around 
binaries with modest eccentricity ^ ~ 0.1 (Artymowicz & Lubow 1994). In Figure 2, we plot the ratio as a function of eccentricity 

in the two cases of a retrograde (T- 1 ^ 2 ) and prograde (ri, 2 ) disc. Notice that the retrograde torque is quite weak compared to the prograde 
torque and there is a strong dependence on eccentricity that is expected to be oc for ^ 1, as seen in the figure. In Figure 3 , we plot 

the torque r_p 2 and the viscous torque at the resonance for the disc parameters adopted in the simulations of Section 4. The dashed line 
shows that r_i 2 follows the expected dependence on eccentricity for small e. For fairly high values fore ~ 0.5, the resonant torque can 
overpower the viscous torque and truncate the disc. 

In summary, of the disc mode strengths we consider, we expect 5 14 , 5 14 and S to be small compared to 1 and 2 . Whereas 
S'-ij can be significant for intermediate eccentricities that are less than about 0.55, it vanishes for equal mass binaries. 5 '_i ,2 should dominate 
among these cases at higher eccentricities ^ > 0 . 6 , especially forq= 1 . 


4 SIMULATIONS 

We present numerical simulations using the smoothed particle hydrodynamics (SPH; e.g. Price 2012) code phantom (e.g. Price & Federrath 
2010; Lodato & Price 2010). This code has already been applied to many types of accretion in binary systems (Nixon 2012; Rosotti et al. 
2012; Nixon et al. 2013; Facchini et al. 2013; Martin et al. 2014a,b). The simulations include the effects of disc pressure and viscosity, while 
effects of disc self-gravity are ignored. The binary is modelled as two Newtonian point masses which affect the gas orbits through gravity 
and accrete any particles which come close enough. The back reaction on the binary orbit is also included, but the disc mass is too small to 
result in any significant changes on the simulation timescale. 


4.1 Setup 

We set up a binary with mass ratio q = M2 1 Mi and total mass M = Mi + M2 - I, semimajor axis a = 1, and eccentricity e initially at 
apocentre of the orbit. The binary is taken to rotate with positive angular speed and so rotate with positive mean motion Qb (defined as 2;r/Pb 
for binary orbital period Pb), while the disc rotates in the opposite sense with D < 0. The binary sink particles have accretion radii of 0.2a, 
inside which gas is removed from the simulation and its mass and momentum added to the sink particle. We initialise the disc with 8 million 
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Figure 2. Torque ratio of the resonant torques T- 1 ^ 2 IT 1^2 at the +1:3 eccentric outer Lindblad resonance for an equal mass binary, assuming the same disc 
density at the resonance in the two cases of a retrograde and prograde disc, respectively. The solid line is the result of numerical evaluation of Equation 21. 
The dashed line follows oc that is expected to be valid fox e 1. 




3 



Figure 3. Upper solid line plots torque T-i^ due to eccentric outer -1:3 retrograde Lindblad resonance in units of at the resonance as a function of 

binary eccentricity e for an equal mass binary. The lower solid line is for a binary mass ratio of 0.1. The dashed line follows oc that is expected to be valid 
for e ^ 1. The horizontal dotted line plots the viscous torque at the resonance for the disc parameters used in the simulations. 


particles in Keplerian orbit about the binary centre of mass. The disc initially extends from an inner radius of la to an outer radius of 8a with 
a surface density distribution of E(r) = Eo(r/ro)"^^^, where the value of So is set to provide a disc mass of = 10“^M. Throughout the 
simulations the binary aceretes < 10% of the disc. The simulations adopt the locally isothermal assumption in which the disc sound speed 
is c^{R) - for spherical radius R and Cs,o corresponds to aspect ratio H/R - 0.05 (0.035) at the disc inner (outer) edge at all 

times. Since the discs are thin, the differenee between between the cylindrical radius r and spherical radius R at any point in the main body of 
the disc is small. We employ an explicit accretion disc viscosity which corresponds to an approximately constant Shakura & Sunyaev (1973) 
a - 0.05 throughout the initial disc (Section 3.2.3 of Lodato & Price 2010). The viscous stresses include a nonlinear term with a coefficient 
/5 av = 2 (AV stands for artificial viscosity) that suppresses inter-particle penetration, as is standard in SPH codes. For this configuration, the 
shell averaged smoothing length per disc scale-height (see Lodato & Price 2010) is (h) IH ^ 0.16 for 8 million particles. The simulations 
reported below span q ^ 0.1,0.3, and 1.0 (more precisely 9:1, 3:1 and 1:1) and e = 0.0, 0.2, 0.4, 0.6, and 0.8. 
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Figure 4. Surface density rendering of each simulation with <5^ = 0.1 at t = 628 (after 100 binary orbits). The binary rotates in a counter-clockwise sense and 
the disc rotates in a clockwise sense. The binary eccentricity is shown on each panel. The colour scheme, from lowest surface density (black) to highest surface 
density (white) covers approximately 2 orders of magnitude. The binary is represented by the two red filled-circles, with the circle size denoting the accretion 
radius inside which particles are removed from the simulation. 


4.2 Results 

In Figures 4, 5, and 6 we show the surface density renderings of the simulations at their completion after 100 binary orbits. For each mass 
ratio the general results are the same: for a circular binary there are no disc structures that signify resonances between the disc and the binary, 
when the binary is eccentric the disc structures appear and generally grow stronger with increasing eccentricity. The unequal mass ratio 
simulations, q - 0.1,0.3, show significantly asymmetric structures, whereas the equal mass ratio simulations (which begin with symmetric 
initial conditions) remain bisymmetric throughout the duration of the simulation. 

The form of the disc disturbances seen in the figures is of spiral waves that have;^^ = r&O^I&r > 0, where 6^{r) is the azimuthal angle of 
the spiral. This sign of x is expected from the theory of disc resonances in Section 3 because each resonance is due to a retrograde bar that 
excites “trailing” waves in the retrograde disc. We note that if the disc were responding to a prograde disturbance, the sign of the;^ would be 
negative. Only if the retrograde disc is responding to a disturbance that rotates in a retrograde sense at a faster speed than the gas can be 
positive, as we find. That condition occurs at an outer Lindblad resonance due to a retrograde bar. In addition, the gas streams in the gap are 
also expected to have x > ^ because the streams approximately conserve angular momentum as they fiow towards the binary. The gas then 
rotates faster in the retrograde sense, resulting inx > 0- 

The circular equal mass ratio simulation displays weak streams feeding the sink particles from the inner edge of the disc. Such features 
may be present in the unequal mass ratio cases, but here the secondary sink is closer to the disc inner edge and accretes the gas which might 
form them. In Figures 7, 8, and 9 we show the mode strengths against time for a selection of the simulations. As the disc is not set up 
in equilibrium with the binary orbit, there is some initial settling of the disc which appears in the mode strengths. Figure 7 shows the five 
calculated mode strengths for the ^ = 0.1, ^ = 0.4 simulation. Figure 8 shows the five calculated mode strengths for the q = 0.3, e = 0.6 
simulation. Figure 9 shows the five calculated mode strengths for the ^ = 1.0, e = 0.8 simulation. These figures illustrate the general trends 
found in this work: (1) smaller mass ratios create weaker resonances, (2) higher eccentricities create stronger resonances, and (3) for strong 
enough resonances the (-1,2) mode can prevent gas from reaching the (-1,1) resonance on circular orbits, depleting its strength. In addition, 
the (-1,1) resonance is likely unstable at higher eccentricties, as discussed in Section 3. 
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Figure 5. Same as Figure 4, but for ^ = 0.3. 


The extent of the gap and the (-1, 2) mode strengths are fairly stable at the end of the simulations. However, since the simulations are 
limited to 100 binary orbits, we cannot be certain of the behaviour over the longer term evolution of the disc. The (l,m) = (-1,2) mode 
strength is the largest and displays the clearest trend of increasing strength with increasing eccentricity, as expected by the analytic model 
in Section 3. The (1,1) mode strength is small for all parameters, while the (1,2) mode strength displays some fluctuations with increasing 
eccentricity. The (-1,1) mode strength nearly vanishes for an equal mass binary and is weaker than the (-1,2) mode strength. The behavior 
of the mode strengths is then in good agreement with the expectations of the analytic model discussed in Section 3. 


4.3 Binary evolution and accretion 

Nixon et al. (2011a) developed an analytic model for the interaction of a binary with a circumbinary retrograde disc. The secondary object 
is taken to be of small mass compared to the primary. The model assumed that the secondary object (or its surrounding disc) experiences 
an impact and accretion from the circumbinary gas that is locally counter-rotating at circular Keplerian speeds. In the model, the binary 
loses angular momentum by gravitational interaction with the retrograde gas. Its semi-major axis decreases, while it eccentricity increases, 
provided that its eccentricity is greater than a threshold value ~ Hjr. 

In the case that the binary opens a gap in a retrograde circumbinary disc, the situation is somewhat different. The binary exerts a tidal 
torque at a resonance, such as the (-1,2) resonance, due to a retrograde component of its potential that is prograde with the rotation of 
the disc, as discussed in Section 3. The resonance has the effect of pushing gas outward, away from binary, which adds negative angular 
momentum to the 12 < 0 disc and therefore adds positive angular momentum to the binary. The resonant torque on the disc is then negative 
and the torque on the binary is positive. In addition, the gas that impacts the secondary flows in the form of gas streams that originate at the 
disc inner edge. These streams are retrograde to the binary and should reduce the angular momentum on the binary. However, they are not in 
circular motion. Therefore, it is unclear how well the model of Nixon et al. (2011a) applies to the case of a highly eccentric binary. 

Here we report the evolution of the binary eccentricity and the mass flow rates on to the sink particles in the simulations. In Figure 10 
we plot the semi-major axis and eccentricity evolution of the binary orbit with time for all of the simulations. In all plots the semi-major axis 
decreases with time, as expected from the capture of angular momentum and energy of retrograde gas orbits. The eccentricity evolution of 
the binary is consistent with the analytical predictions of Nixon et al. (2011a). The circular binaries (for which e ^ Hjr) remain circular as 
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Figure 6. Same as Figure 4, but for q = 1.0. 



Figure 7. The different mode strengths as a function of time for the q = 0.1, e = 0.4 simulation. The time axis is in units where Itt is one binary orbital period. 
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Figure 8. The different mode strengths as a function of time for the q = 0.3, e = 0.6 simulation. The time axis is in units where In is one binary orbital period. 



0 200 400 600 

time 


Figure 9. The different mode strengths as a function of time for the < 5 ^ = 1.0, e = 0.8 simulation. The time axis is in units where 2n is one binary orbital period. 


accretion at apocentre of the binary, which increases eccentricity, is offset by accretion at pericentre, which decreases eccentricity. For these 
simulations the eccentricity varies little and is approximately 10""^. The model of Nixon et al. (2011a) predicts that the eccentricity e should 
grow if it satisfies e> Hjr, and our simulations generally confirm this. Some of the simulations show eccentricity decay at early times before 
the disc has settled into a quasi-steady state. The only simulation which does not follow the predicted trend is ^ = 1.0, ^ = 0.8. For this case 
the eccentricity decays from 0.8 to 0.799 over the timescale of the simulation. It is unclear if the eccentricity would continue to decrease if 
the simulation were run for longer, or instead later increase. 

The remaining eccentric simulations all show eccentricity growth with dejdt ~ 10“^ in units of time where 2n is one binary orbital 
period. This is in agreement with the analytical prediction of ~ M/M 2 of Nixon et al. (2011a) (see Figures 11, 12, and 13 for the accretion 
rates). 
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Figure 10. Evolution of the binary semi-major axis, a, and eccentricity, e. The time is in units where In is one binary orbital period. The semi-major axis is 
shown on the left hand panel and the eccentricity is shown in the right hand panel. From top to bottom the rows show the different initial eccentricities from 
^ = 0.0 - 0.8. The time resolution of the plot is one data point every 10 binary orbits. 
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Figure 11. Accretion rate of gas flowing on to the binary for the ^ = 0.1 simulations. This measures the mass flux through the sink particle radii. The accretion 
is binned with width equal to one binary dynamical time, 1/Qb- The accretion is the total amount on to the primary and secondary sinks. Only the last 10 orbits 
of the binary are shown for clarity. 



In Figures 11, 12, and 13, we plot the accretion rate of gas on to the two sink particles as a function of time. The curves are binned with 
width equal to one binary dynamical time, l/Hb. We do not separate the primary and secondary accretion rates, as the secondary accretes 
almost all the gas, except for the equal mass case, where the rates are similar. In general the accretion rate is similar to the circular binary 
case, but with an obvious binary orbital period modulation. This modulation has an amplitude of approximately 1.5 orders of magnitude for 
^ = 0.1 and ^ = 0.3, but up to 3 - 4 orders of magnitude for the equal mass case. This periodicity is marked by strong accretion at apocentre 
and little accretion at pericentre of the binary orbit. In the prograde case, there is also a strong periodic modulation of accretion rate over the 
binary period for eccentric orbit binaries. However, the accretion rate is highest near periastron (Artymowicz & Lubow 1996). 


5 DISCUSSION 

We first discuss the limitations of our simulations, and possible improvements. We then consider the consequences of our results for possible 
astrophysical applications. 

The simulations described above employ a simple locally isothermal thermodynamic treatment for the gas, where the sound speed is a 
radial power-law from the centre of mass. This approximation should not alter the fundamental behaviour of the discs, but does prevent any 
thermodynamic effects from arising. It would be interesting to run these simulations with an equation of state that accounts for shock and 
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viscous heating while allowing the gas to cool on an appropriate timescale. This may affect the propagation of waves (e.g. Lubow & Pringle 
1993) and thus the redistribution of angular momentum by the resonances. 

Similarly our simulations neglect the effects of self-gravity in the gas, and the effects of magneto-hydrodynamic turbulence induced by 
the magneto-rotational instability, instead employing an accretion disc viscosity through an a parameter. These additional physical processes 
should be included in a full treatment of the problem, but are not expected to significantly alter our conclusions. 

The simulations could also be run for longer to discover longer timescale evolution of the disc structures. However, the limiting factor 
is that accretion and disc spreading can deplete the gas to the point where there is too little of it in regions of interest. Therefore it would be 
worthwhile including mass input into the disc to allow a quasi-steady state to be achieved over a long timescale simulation. 

Binaries exist in a variety of astrophysical scenarios. However, to create a retrograde planar disc probably requires a chaotic environment 
which can introduce gas to the binary with random orientations (but see also e.g. Pringle 1996, 1997 for a tilting instability which can cause 
the tilt of an initially planar disc to grow to retrograde angles). The two most likely chaotic scenarios are stellar binaries formed in dense 
star forming regions (e.g. Bate et al. 2010) and SMBH binaries accreting from the host galaxy (e.g. King & Pringle 2006, 2007; King et al. 
2008). In both cases, a significant fraction of randomly oriented discs are expected to align or counteralign with the binary plane depending 
primarily on the initial disc-binary inclination angle, but also on the ratio of disc to binary angular momentum (King et al. 2005; Nixon et al. 
2011b). 

Retrograde accretion from circumbinary discs has recently been proposed by Nixon et al. (2011a) as a possible solution to the last parsec 
problem (Begelman et al. 1980; Milosavljevic & Merritt 2001). This model relied on the eccentricity growth, promoted by the capture of 
negative angular momentum gas, continuing to drive the binary eccentricities to large values, forcing the binary into the gravitational wave 
regime. If the resonances discussed here for retrograde circumbinary discs inhibit eccentricity growth as is expected in the prograde case 
for large eccentricities, then this may not be possible. However, these resonances require significant eccentricity to achieve large strengths, 
and even so are of significantly smaller strength than their prograde counterparts. Therefore it is quite possible that efficient eccentricity 
growth to values approaching unity is still possible. For this discussion the accretion rates on to the binary are highly suggestive, as the mean 
accretion rate appears to be insensitive to the binary eccentricity (see Figs 11, 12, and 13). Instead the eccentricity acts to introduce significant 
periodicity in the mass flow rates through the sink accretion radii. As the binary is still capturing the same amount of material each orbit, it 
appears likely that the conclusions of Nixon et al. (201 la) are independent of the presence of these resonances. This also appears likely from 
Figure 10 as the eccentricity generally grows in our simulations and at a rate consistent with the analytical prediction of Nixon et al. (201 la). 
However, as discussed in Section 4.3, the assumptions of the Nixon et al. (2011a) model may not be well satisfied at high eccentricity. We 
shall revisit this issue in more detail in future work. 


6 CONCLUSIONS 

We have presented an analytical model and three-dimensional hydrodynamical simulations of retrograde circumbinary discs covering a 
range in binary mass ratio and eccentricity. These simulations are consistent with the analytical prediction that binary eccentricity causes 
disc resonances to occur, which leads to angular momentum exchange between the binary and the disc. For circularly orbiting binaries, there 
are no disturbances (modes) excited at Lindblad resonances in a retrograde disc. But for eccentric orbit binaries, Lindblad resonances can 
be excited in such a disc. The Lindblad torques increase with eccentricity. At high eccentricity ^ ^ 0.6, the most strongly excited disc mode 
(-1,2) found in the simulations is also the mode expected to be most strongly excited in the analytic model by the Lindblad resonance 
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at Q(r) ^ -^b/S. The Lindblad torques associated with this mode can open a gap in the disc at binary eccentricity e^O.6 for typical 
disc parameters. The inner resonances are destabilized at high eccentricities. We also find that the binary semi-major axis and eccentricity 
evolution are in general agreement with the analytical predictions of Nixon et al. (2011a). 

The binary eccentricity evolution is affected by gas accretion and resonant interactions with the disc. Future simulations will focus on 
whether the binary eccentricity can grow to approximately unity or if there is a smaller limiting value. By including mass injection into the 
disc it will be possible to evolve the discs to a quasi-steady state and explore longer term evolution than is described here. 
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Figure Al. Surface density rendering of the same simulation after 100 binary orbits, shown at three different resolutions. The colour scheme represents 2 
orders of magnitude in surface density from the highest density (white) to the lowest density (black). The binary is represented by the two red filled circles, 
where their size denotes the size of their accretion radii. The three pictures are shown at precisely the same run time (100 binary orbits). The different phase of 
the binary orbit is due to the slightly different back reaction from the disc in each case. 



Figure A2. Comparison between low, medium and top resolution simulations for q = 0.3. The black-solid line corresponds to the lowest resolution (1.25 x 10^ 
particles), the red-dashed line corresponds to the medium resolution (10^ particles) and the green-long-dashed line corresponds to the highest resolution 
(8 X 10^ particles). For e = 0 the mode strength should be zero and therefore the magnitude drops with increasing resolution. As the eccentricity is increased 
the mode strength increases. For the medium and highest resolution simulations good agreement is found, but at the lowest resolution it is clear the modes are 
not sufficiently resolved. 


APPENDIX A: RESOLUTION 

Each simulation was run at three different resolutions corresponding to 1.25 X 10^, 10^ and 8 x 10^ particles. At each different resolution, 
the SPH viscosity was scaled so as to give the same physical viscosity ass = 0-05 using the method described in Section 3.2.3 of Lodato & 
Price (2010). In Figure Al we show the surface density rendering of the ^ = 0.3, ^ = 0.6 simulation at each resolution. It is clear from this 
figure that the lowest resolution simulation is under-resolved, but that at higher resolution the behaviour shows convergence. In Figure A2 
we show the (-1,2) mode strength for each eccentricity at the different resolutions. This shows the same picture, where the lowest resolution 
simulations are not resolved, but the medium and highest resolution simulations show good agreement. For ^ = 0 the mode strength drops by 
approximately the same factor as the resolution is increased - this is to be expected as the mode strength should be formally zero. For higher 
eccentricities the difference between the medium and high resolution simulations is small. 
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APPENDIX B: APPLICATION OF MODE STRENGTHS TO SPH SIMULATIONS 


In the application to SPH simulations, we consider N particles each of unit mass that are located at the 3D Cartesian positions (x/r), ^/r), Zjir)) 
and approximate the surface density distribution as 

N 

2(x, y,T) = Yj - Xj(T))6(y - yj(T)) (B1) 

7=1 


for a disc of mass N. In general, we have 

= I ^f,g,e,mir)rdr 
Jo 
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= 0 2 n X I I I ^(r,d,T)f{m0)g{£T) rdrdddT 
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f(m6j(T))g(£T)dT, 


where / and g each can be the cos or sin functions and Ojir) = arctan (yj(T)/xj(T)). 

As a illustration of this analysis, we consider a set of N particles labelled by j at radius r = 1 having angular location 
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where 
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for integer 1 < j < N and |/I| «: 1. It then follows that the density distribution of particles is 
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Eor ^ ^ 0 or m ^ 0, we have 
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where the -i factor is included to obtain the sin(mi^y - £t) function in Equation (BIO). 

We show below that this expression for can be also obtained from the use of Equations (6) and (B5). We use the fact that 


sin {mOj) = sin {mil/j + mX cos (mif/j - £t)), 
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since \A\ ^ 1. We then evaluate 
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(B20) 
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From equation (6), we have 
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that agrees with Equation (Bll) and so analytically shows the validity of the decomposition prescription. We also have that 
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